ACID RAIN: MESOSCALE MODEL 


By Hsiao-ming Hsu* 


INTRODUCTION 


A mesoscale numerical model of the Florida peninsula has been formulated 
and applied to a dry, neutral atmosphere. The preliminary results of this 
model presented in this paper are very encouraging. To simulate synoptically 
undisturbed and disturbed situations in this area, several modifications have 
appeared desirable; thus, a new model quite different from the one just men- 
tioned has evolved. A submesoscale model for the Cape Canaveral area will 
use the results from the mesoscale model as continuous inputs. The prospec- 
tive use of the Control Data Corporation (CDC) STAR-100 computer at the NASA 
Langley Research Center for the submesoscale model is discussed. 


PRELIMINARY RESULTS 


The numerical model presented in reference 1 is tested under synoptically 
undisturbed conditions. Because this test is for demonstrating the simulation 
of only the basic dynamics in the mesoscale region and for detecting any pro- 
graming error in the model, three assumptions are applied to simplify the 
calculation. 

1. The atmosphere is assumed to be dry; i.e., the parameterization of 
cloud and precipitation physics is not yet included. 

2. A constant potential temperature (0 = 298.16 K) and a constant 
geostrophic windspeed (V = 6 m/sec) are assigned to the model as the initial 
conditions . 

3. The daytime land surface temperature is represented by a single sine 
function that limits the maximum heating to 10 K at the end of the sixth hour 
of simulation. The water surface is assigned a constant temperature that is 
the same as the initial temperature of the land surface. (The time step is 

5 minutes.) 

Two cases, differing only in the direction of the prevailing geostrophic 
wind, are examined: a prevailing southwest wind and a prevailing southeast 

wind, both 6 m/ sec at all levels initially. 
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Prevailing Southwest Wind 


Results from the simulation experiment of prevailing southwest wind are 
shown for the various levels and times in figures 1 to 3. The horizontal 
wind field at the 50~meter level shows only a slight change from the prevail- 
ing wind at the end of 2 hours (figs. 1(a) and 1(b)). By the third hour 
(fig. 1(c)), the low-level wind direction has turned noticeably at the shore- 
lines and especially sharply along the east coast of the Florida peninsula. 

A definite low-level convergence zone is formed by the end of the fourth hour 
of simulation (fig. 1(d)). The sea-breeze convergence zone continues to de- 
velop and moves inland during the fifth (fig. 1(e)) and sixth (fig. 1(f)) 
hours . 

The horizontal wind fields at 200-, 450-, 800-, and 1250-meter levels are 
shown at the end of the sixth hour in figure 2. The change of wind from 50 
to 1250 meters reveals the vertical sea-breeze circulation along the east 
coast: the sea breeze gives onshore winds in the lower levels and offshore 

winds in the higher levels. The gradual veering of the wind from southeast 
or east near the surface to southwest aloft is clearly shown. 

Figure 3 gives the computed horizontal distributions of vertical velocity 
at 1250 meters at the end of each hour of simulation. The velocity fields 
are mainly developed by the effects of advection over the gradually warmed 
peninsula. The convergence zone forms first along the west coast, whereas di- 
vergence prevails along the east coast. By the third hour, convergence is 
also noted along the east coast. Several distinct maximums of positive verti- 
cal velocity exist along the east coast at the sixth model hour. The maximum 
south of Lake Okeechobee, predicted by Pielke (ref. 2), has the greatest upward 
velocity. Two other maximums north of Lake Okeechobee, although weaker than 
the one just mentioned, are very close to the area of Cape Canaveral. Pielke 
has demonstrated that the predicted regions of positive vertical velocity 
agree qualitatively with radar and satellite observations of cloud and shower 
activity over southern Florida on undisturbed days. His study did not, how- 
ever, extend north of Lake Okeechobee and so did not reveal the effects at 
Cape Canaveral. The maximums of positive upward motion in the author's re- 
sults are expected to be closely related to cloud and shower activity, at 
least qualitatively. In this case, having the southwest wind in a synoptical- 
ly undisturbed situation, disturbed weather contributed solely by the sea- 
breeze convergence appears both in the Cape Canaveral area and in southeastern 
Florida. 


Prevailing Southeast Wind 

The horizontal wind patterns for prevailing southeast wind at the 50- 
meter level are shown in figure 4 for the end of each model hour. Again, a 
slight turning is seen at the coastlines by the end of the second hour (fig. 
4(b)). The convergence zone develops gradually along the west coast of 
Florida. At the sixth hour, the horizontal wind distributions at heights of 
200, 450, 800, and 1250 meters indicate the sea-breeze circulation along the 
west coast (fig. 5). The vertical velocity fields at the 1250-meter level 
at the end of each model hour (fig. 6) show the evolution of the sea-breeze 


34 



convergence zone. Under the southeast wind regime, the area of Cape Canaveral 
is covered by a field of downward vertical velocity after the second hour 
(fig. 6(b)), but a band of upward motion that occurs inland during the fifth 
and sixth hours (figs. 6(e) and 6(f), respectively) indicates possible weather 
disturbance. Along the western shore of Florida, however, the vertical veloc- 
ity fields are more strongly developed. 

For both cases, the motion fields appear to be somewhat exaggerated 
after 5 hours of simulation. This result is a product of the highly idealized 
assumption that the atmosphere is adiabatic (constant potential temperature 
0 = 298.16 K) initially. The computations are stopped at the end of the 
sixth hour. 


Discussion 

Analogous features are shown under the two different prevailing wind 
directions. They are not caused by different dynamics but by the orientation 
of the Florida peninsula with respect to the wind. On synoptically undis- 
turbed days, the local weather is dominated by the sea breeze. For the Cape 
Canaveral area, locally disturbed weather can be expected to appear under a 
general southwest flow. Under a general southeast flow, disturbances at the 
Cape are not expected, but they might develop inland and downwind from the 
Cape. This conclusion agrees qualitatively with Neumann's report that the 
probability of having afternoon thunderstorms under southeast winds is low 
but that the probability of having them under southwest winds in the summer 
months at the Cape is more than 50 percent (ref. 3). 

Comparison of the author’s results with those of Pielke (ref. 2) re- 
veals that the strength of upward motion is weaker in the author's simulation 
than in Pielke 's before the end of the fifth hour. Because the author's hori- 
zontal grid spacing is about three times longer than Pielke 1 s, the extremes 
shown by Pielke' s model are necessarily greater than those shown by the 
author's. At the same time, the author's model is flexible in that the grid 
spacing can be adjusted in both the horizontal and vertical directions. The 
question of increased resolution of the model in turn raises the question of 
the limits of computer capacity. The author proposes, as a second stage of 
model development, the introduction of a cloud-scale model (grid spacing 2 
by 2 by 0.5 kilometers) the initial and boundary conditions of which will be 
derived from the mesoscale model. 

Although it is clear that the modeled mesoscale dynamic processes do pre- 
dict the sea-breeze regimes along the Florida coasts, it is also clear that 
important modifications must be introduced to achieve more complete simula- 
tions. These are presented and discussed in the following section. 
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RECENT DEVELOPMENTS 


To generate more complete mesoscale information over the region of the 
Florida peninsula, the following additional considerations are introduced. 


Parameterization of Convective Cloud and Precipitation 

The study of parameterization described in the section entitled "Future 
Work" leads to the conclusion that Kuo f s method (ref. 4) is the best approach 
for meeting the requirements of fidelity with simplicity. This parameteriza- 
tion is recommended partly by its more rigorous derivation and partly by its 
potential for further development. It will, for example, enable the use of a 
convective cloud model to obtain values for the cloud temperature and specific 
humidity. 

The parameterization allows for two kinds of moisture sources for a con- 
vective cloud (fig. 7). One is the net convergence of moisture between cloud 
base and cloud top 


~f C [fc Cp<lu> + |? (pqv) +|^Cpqw)]d Z 
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and the other is the evaporation between the surface and cloud base 
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where z t and z^ are the heights of cloud top and cloud base, respec- 
tively; u and v are the x- and y-component velocities; w is verti- 
cal velocity; q is the specific humidity; p is density; Cp is the drag 
coefficient; V is the horizontal windspeed; and g denotes the surface. 

The total moisture available M may fuel the convective cloud by becoming 
entrained through the cloud sides, by pumping moisture through the cloud base, 


36 



or by evaporating part of M in the environment. Two kinds of sink from a 
convective cloud are parameterized: the evaporation from cloud sides and 

the precipitation received at the surface. The evaporation both in the en- 
vironment and from cloud sides represents moisture input to the air column 
8M. The precipitation received at the surface represents the removal of 
moisture (1 - 0)M from the air column. The processes within the dashed box 
in figure 7 are the subgrid processes in the author's mesoscale model and 
are not evaluated explicitly. 

The source term for the potential temperature 0 can be expressed by 



( 2 ) 


where L is the latent heat of evaporation and C is the precipitation rate 
produced by a convective cloud. The value it is the scaled pressure 


TT 


R 




(3) 


where Cp is the specific heat at constant pressure, R is the gas constant 
for dry air, p is the pressure, and p 0 = 1000 millibars is a reference 
pressure. 

The vertically integrated precipitation rate gives the total amount of 
precipitation. 



pc dz = (1 - 0)M 


(4) 


The precipitation rate produced by a convective cloud may be written as 


pC = (1 - g)M 


N (z) 


K - %) 


(5) 
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if the vertical distribution function N(z) of 0 satisfies 



( 6 ) 


Supposing that N(z) is represented as a function of the difference of po- 
tential temperatures between the cloud interior and its environment, 


N(z) 


( 9 c - 9 ) 

<e - e> 

c 


(7) 


where 


<0 

c 



( 8 ) 


Then, 


(1 - p)M ( 6 c 9 ) 
p ( Z t " Z b) <9 c ‘ 0> 


(9) 


where the subscript c denotes cloud. 

Substituting equation (9) into equation (2) gives 

= L(1 - B)M ( 9 c ~ 9 ) 

b 02 TTp^Z t - <0 c - 0> (10) 
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for the source term for 0. 


Similarly, the source term for q is 


s = (1 - 8)M ( e c e ) 

q2 p ( z t " z b) <6 c - 0> 


( 11 ) 


The coefficient 3 representing the fraction of M remaining in the air 
column is a very important and crucial factor. The whole evaporation process 
is controlled by 3- It has been assumed to have the following form (ref. 5). 


- (i/q, 


-.n 


3 = 


L 1 ' ( q / q s )*J 


when (q/q s ) > (q/q g )* 


1, when (q/q s) < (q/q s) * 


( 12 ) 


where the subscript * represents a critical value to be specified and the 
superscript n is a constant. 

The cloud variables that have the subscript c are determined by the 
saturated adiabat between cloud base and cloud top. It is intended that 
these values will eventually be computed by using a convective cloud model. 


Diagnostic Vertical Velocity Equation 

It is intended to represent all diurnal, seasonal, and synoptic character- 
istics of the region of interest. A shallow model clearly does not meet this 
requirement. To account for deep convection, the model is extended to 16 kil- 
ometers vertically in 18 levels (table I). The assumption of constant den- 
sity can no longer be used and must be replaced by the complete continuity 
equation for a compressible fluid. Finally, a modified Richardson's equation 
is used to calculate the vertical velocity diagnostically. 

The principle of mass conservation gives the continuity equation in 
three dimensions. 


V*V 


p dt 


(13) 
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where 


u ^ t , 3 *t , 9 ^ 
V = — l + — 1 + — k 
3x 8y J 8z 


V = 


ui + vj + wk 


and 


Taking the logarithm of the equation of state and the definitions of 
IT and 0 yields 


to p = Jin p + to R + to T (14) 

to 7T = to c + R/c top - R/e to p n (15) 

p / p / p u 

Jin 0 = k T - £n tt + £n c (16) 

P 


where T is the temperature. Differentiating equations (14), (15), and 
(16) and making substitutions, one obtains 


1 _ jdp _ Sr 1 _ djr 1_ (10 

p dt R it dt 0 dt 


where c v is the specific heat at constant volume. Substitution of equation 
(17) into equation (13) gives 


V*V 


C v _1 tin | 1 d6 
R it dt 0 dt 


(18) 
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Multiplying p/ff on each side of equation (18), taking ^*Vinr in dir/dt 
from the right-hand side of equation (18) to the left-hand side, and 
rearranging terms, one has 
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R 2 9t 
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The time derivative of the hydrostatic equation 
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is integrated from any level z to the model top H to determine 
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Both 3/3t(l/0 v ) and d9/dt can be evaluated explicitly from the equations 
of potential temperature and specific humidity in the model. Then, equation 
(19) becomes 


V* 



“ ~ ^2 (A + B) + 

TT 


TT© dt 


(24) 


after substituting equation (21) into equation (19). 

Integrating equation (24) from the surface to any height z and 
applying the surface boundary condition for w, w * 0 at z = 0, one 
has 


“ ■ ■ p/ [fj(» “ ) + ^ V )] dz 


/»2 

-- f 
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< A + B > -fez? 

TT 


dz 


(25) 


This equation gives the values of the vertical velocity at any point and 
is called the "modified" Richardson's equation. 
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If the top boundary condition for w, w = 0 at z = H, is further 
applied, the pressure tendency will be given as follows. 



(26) 


The pressure is adjusted by the predicted value from equation (26) at the 
model top, and the rest of the TT-field is calculated from the hydrostatic 
equation. 


Lateral and Upper Boundary Conditions 

The Neumann-type lateral boundary condition gives some noise along the 
boundaries. It can be found in figures 3 and 6. The amplitude of this noise 
is small, but the noise will be amplified in a stably stratified atmosphere 
and become fast-moving internal gravity waves. These waves can seriously 
contaminate the model results. Perkey and Kreitzberg (ref. 6) developed a 
method to control this problem; this method has been adapted to the present 
model. 

For both lateral and upper boundaries, the variables u, v, 9, 
and q (represented by D) on the boundary region will be determined by 
weighting their local tendencies. 
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where the subscripts m and % represent mesoscale and large scale, 
respectively. The weighting functions for 3D/3t are as follows. 


W 1 (B) 

W (B ± 1) 
W 1 (B ± 2) 
(B ± 3) 
W 1 (B ± 4) 



( 28 ) 


where B denotes the grid point at the boundary, B ± 1 is the grid point 
next to B (one grid space in from the boundary), etc. 

The vertical velocity w is also determined by means of a weighting 
procedure as follows. 


W(I) = W 0 (I) w(I) 

z m 


(29) 


The weighting functions for w are as follows. 


W 2 (B) 

w 2 (b ± 1) 

W 2 (B ± 2) 
W 2 (B ± 3) 
W 2 (B ± 4) 



(30) 


These weighting functions may be adjusted, depending upon test results. 
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Boundary-Layer Parameterization 

To understand the full diurnal variation of mesoscale weather is one of 
the author's goals. Deardorff's equation (ref. 7) may be used to model the 
daytime growth of the height of the planetary boundary layer (PBL). The equa- 
tion is not, however, capable of simulating the variation of the height of 
the PBL after sunset and before sunrise next morning. To predict the height 
of the PBL at night, the formulation developed by Smeda (ref. 8) will be 
used. 


8H 
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( 9H * 
Tk i^ + v 


H * 3 y 


+ 0.06 


H * £ 


1 - 
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(31) 


where f = 2ft sin <J ), ft is the Earth rotation rate, <)) is the latitude, 
is the height of the PBL, and u. is the frictional velocity. The vertical 
profile of the eddy exchange coefficients remains the same as before. 


FUTURE WORK 


Improvements of the Mesoscale Model 

The mesoscale model is now quite well developed and will be used to pro- 
vide basic predictive information for the Florida peninsula. As is evident 
in the model derivations, it has been necessary to idealize and/or parameter- 
ize numerous component parts of the problem. It is anticipated that develop- 
ment of techniques and procedures for handling these component parts will 
continue and that improvements of the model will be made possible by these 
developments. This is a continuing process which will provide improvement 
but cannot be expected to achieve perfection in the model. 


Development of a Submesoscale Model 

To simulate the processes of the solid rocket motor exhaust product re- 
moval and the patterns of deposition that they produce in the area of Cape 
Canaveral after Space Shuttle launches, a submesoscale numerical model is being 
developed. The model domain is in three dimensions, 90 by 90 by 15 kilometers. 
The approximate location is shown in figure 8. Processes will be depicted 
in terms of the behavior of the meteorological variables and model parameters 
at some 65 596 grid points (a spacing of 2 by 2 by 0.5 kilometers). The basic 
equations are derived from the theory of atmospheric dynamics and thermodynam- 
ics, cloud and precipitation physics, boundary-layer dynamics, and air pollution 
chemistry. Initial and boundary conditions will be obtained continuously 
from the output of the mesoscale model. 
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The requirements of this model in terms of computer capacity and speed 
are very great. Consequently, a comparative study of the available and soon- 
to-be available computer facilities has been undertaken. 


COMPARATIVE STUDY OF COMPUTER FACILITIES 


A brief comparative study of the prominent advanced computer systems 
that are or may become available has been made. The systems considered are 
the CDC-6600, the CDC-7600, the Amdahl-470V/6 , the CRAY-1, and the STAR-100. 

Central processing unit (CPU) rates and memory storage capacities are 
tabulated and compared in terms of performance ratios in table II. The internal 
core memory is augmented by virtual memory of unlimited capacity in the Amdahl 
and the STAR systems. This may be a necessary feature for the submesoscale 
model. 

Clearly, the STAR-100 system is comparable to the CRAY-1 in terms of 
the CPU rates. The variation of these rates for a given unit is a function 
of the kind and extent of the computations and tends upward for the array com- 
puters as the number of grid points increases. For the STAR-100, in particu- 
lar, the CPU performance ratio with respect to the CDC-6600 varies as shown 
in table III. 

It is, of course, not valid to extrapolate the trend shown. The antici- 
pated model space of 65 596 grid points will clearly require use of the vir- 
tual memory capability of the STAR. At present, no means exists for estimat- 
ing how the size of the model will affect the overall speed of simulation. 

Advantages of the Amdahl-470V/6 at the University of Michigan and of the 
CRAY-1 /CDC-7600 at the National Center for Atmospheric Research (NCAR) are 
that both have sophisticated graphics facilities and that the CRAY has both 
scalar and vector registers. It is not yet clear to what extent these fea- 
tures might simplify the procedures and facilitate the presentation of simu- 
lation results. 

One consideration of significance in selecting the computer facility to 
be used is that of competing demands for computer times. The NCAR is cur- 
rently inviting requests for time on the CRAY, and there is little doubt that 
the atmospheric modeling community will respond. If the STAR is made avail- 
able on comparable financial terms to the present project, the much more 
restricted competition for its use will probably make it the more attractive 
alternative . 
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TABLE I.- VERTICAL GRID FOR DEEP CONVECTION MODEL 
HEIGHT OF VERTICAL LEVEL 


Level no. 

Height, km 

18 

16.00 

17 

14.00 

16 

12.00 

15 

10.00 

14 

8.50 

13 

7.00 

12 

6.00 

11 

5.00 

10 

4.05 

9 

3.20 

8 

2.45 

7 

1.80 

6 

1.25 

5 

.80 

4 

.45 

3 

.20 

2 

.05 

1 

0 
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TABLE II.- COMPARATIVE CPU RATES AND STORAGE CAPACITIES 


a Amdahl-470V/6 b CDC-6600 b CDC STAR-100 c CDC-7600 

c CRAY-l 

CPU performance ratio 



1.0 

to 

1 . 

,33 

0.37 

2.59 

to 24.89 

1.85 

3.7 to 14.8 

2.7 

to 

3. 

,6 

1.0 

7.00 

to 67.2 

5.0 

10.0 

to 

40.0 

0.39 

to 

0 . 

.51 

0.14 

1.00 

to 9.6 

0.71 

1.43 

to 

5.71 

0.04 

to 

0 . 

,05 

0.01 

0.10 

to 1.0 

0.07 

0.11 

to 

0.42 

0.54 

to 

0 , 

.72 

0.2 

1.4 to 13.44 

1.0 

2.0 

to 

8.0 

0.27 

to 

0 , 

.36 

0.1 

0.7 

to 6.72 

0.5 

1.0 

to 

4.0 

0.03 

to 

0 , 

.09 

0.03 

0.175 

to 1.68 

0.13 

0.25 

to 

1.0 


CPU word storage capacity 


1 048 576 131 000 512 000 512 000 1 048 576 

plus VM d plus VM d 


a University of Michigan. 
b NASA Langley Research Center. 
c National Center for Atmospheric Research. 
^Virtual memory. 


TABLE III.- CPU PERFORMANCE RATIOS FOR STAR-100/CDC-6600 
IN RELATION TO NUMBER OF GRID POINTS 


No. of grid 
points 

CPU performance ratio 
STAR-100/CDC-6600 

168 

13.5 

280 

12.3 

336 

20.3 

504 

24.4 

2800 

67.2 
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(a) After first hour. 


(b) After second hour. 
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(c) After third hour. 


(d) After fourth hour. 


Figure 1.- Horizontal wind field for various 
Velocity vectors are in meters per second, 
undisturbed, neutral, and dry air; initial 
southwest at all levels.) 


elapsed times at 50-meter level 
(Conditions are as follows: 
windspeed of 6 m/sec from 
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(e) After fifth hour. 


(f) After sixth hour. 


Figure 1.- Concluded 0 
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(a) After first hour. 


(b) After second hour. 


Figure 3.- Distribution of vertical wind velocity (centimeters per second) 
for various elapsed times at 1250-meter level. (Conditions are the same as 
for fig. 1 . ) 
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(c) After third hour. (d) After fourth hour. 

Figure 3.- Continued. 
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(e) After fifth hour. 


(f) After sixth hour 


Figure 3.- Concluded 
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(a) After first hour. 


(b) After second hour. 
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(c) After third hour. 


(d) After fourth hour. 


Figure 4.- Horizontal wind field for various elapsed times at 50-meter level. 
Velocity vectors are in meters per second. (Conditions are as follows: 
undisturbed, neutral, and dry air; initial windspeed of 6 m/sec from 
southeast at all levels.) 
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(e) After fifth hour. 


(f) After sixth hour. 


Figure 4*- Concluded * 
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(a) At 200-meter level* 


(b) At 450-meter level. 
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(c) At 800-meter level. 


(d) At 1250-meter level. 


Figure 5.- Horizontal wind field at various altitudes after sixth hour. Ve- 
locity vectors are in meters per second. (Conditions are the same as for 
fig. 4.) 
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(a) After first hour. (b) After second hour. 

Figure 6.- Distribution of vertical wind velocity (centimeters per second) 
for various elapsed times at 1250-meter level. (Conditions are the same as 
for fig. 4.) 
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Figure 7.- Parameterization scheme for convective cloud and precipitation in 

the mesoscale model. 
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Figure 8.- Map showing the area to be covered by the submesoscale model. 


63 




